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(^ Abstract 

(N 



This paper is concerned with the numerical minimization of energy func- 
^' tionals in BV{Q) (the space of bounded variation functions) involving total 



variation for gray-scale 1-dimensional inpainting problem. Applications are 
shown by finite element method and discontinuous Galerkin method for total 
variation minimization. We include the numerical examples which show the 
different recovery image by these two methods. 
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QN^ 1. Introduction 



In the first chapter of the book \T\ Holger Rauhut has already intro- 
O duced that the minimization of £i-norms occupies a fundamental role for the 

promotion of sparse solutions. This understanding furnishes an important in- 
terpretation of total variation minimization |2j as a regularization technique 
for image inpainting. In this paper we consider as in [3l H] the minimization 
in BV{Q) (the space of bounded variation functions O |6]) of the functional 



J{u):= f \Tu{x)-g{x)fdx + 2\\Du\{n), (1) 

Jn 

where Q C R^, for rf = 1,2 be a bounded Lipschitz domain, T : L^{Q) — )■ 
L'^{Q) is a bounded linear operator , g G L'^{Q) is a datum, \Du\ {fl) := 
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J^ |V'u(x)| dx is the total variation of u, and A > is a fixed regularization 
parameter [2| . Several numerical strategies to efficiently perform total varia- 
tion minimization have been proposed in the literature, refer to [8[ |9], [101 HI] ■ 
However, the interesting solutions may be discontinuous, e.g., along curves 
in 2D. Hence, the crucial difficulty is the correct numerical treatment of 
interfaces, with the preservation of crossing discontinuities and the correct 
matching where the solution is continuous instead, see Section 7.1.1 in [T2] . 
In order to deal promptly with the discontinuity, we have studied the appli- 
cations to gray-scale 1-dimensional inpainting problem by the finite element 
method and discontinuous Galerkin method(Refer to [131 [H]) for total vari- 
ation minimization, respectively. 

The paper is organized as follows. Section 2 is devoted to the alternating- 
minimization algorithm to compute minimizers of J{u). In Section 3, the 
method of finite element method for total variation minimization for our 
problem is illustrated. The main work of this paper about discontinuous 
Galerkin method for total variation minimization is described in Section 4. 
Finally, we include some numerical experiments and discuss their results and 
describe our future study. 

2. Euler-Lagrange equation and a relaxation algorithm 

In this section we propose a method for solving the total variation min- 
imization problem ([I]) in 1-dimensional case. The details could be found 
in [15]. For gray-scale 1-dimensional inpainting problem, the functional ([T]) 
becomes 

J{u):= I \Iq:\d{u{x) - g{x))\ dx + 2\ / \u'{x)\dx, (2) 

Jn Jn 

where D C fl is the damaged domain with measure fi{Q \D) > 0, and 1q\d 
denotes the characteristic function oi Q\D. 

Associated to J^ we have the formal Euler-Lagrange equation: 

u' 
-X( — y + (u~g)ln\D = 0, (3) 

with suitable boundary conditions. In our case, we use Neumann conditions. 
Later we introduce a new functional given by 

eh{u,w)=2 / \1q\£){u{x) — g{x))\ dx + 2X {w\u'\ H )dx, (4) 

Jn Jn w 



where u G W^'"^ {fl; R) , and w G L'^{Q;R) is such that eh < w < — , where 
{eh} is a positive decreasing sequence such that lirrih^ooeh = 0. While the 
variable u again is the function to be reconstructed, we call the variable w 
the gradient weight. 

For any given u^'^^ and w^^\ we define the following iterative alternating- 
minimization algorithm: 



y{n+i) ^ argmmuew^,2^n.R)e{u,w^''^), 

)("+!) = ar5fmin<^< j_ £('u'^"'+-^\t(;). '^ '' 



M7^ 



Then we have the 1-dimensional convergent result of Theorem 7.2 in [15]. 

Theorem 1. The sequence {«'•"■' }„gN /^Qs subsequences that converge strongly 
in L^(r2; M) anc? weakly in W^''^{^; M) to a stationary point u^°°^ of J'; i.e., 



(oo) 



solves the Euler- Lagrange equations M). Moreover, if J has a unique 
minimizer u* , then u^°°^ = u* and the full sequence {^^'■"''jneN converges to 



u^ 



From Theorem IT] we conclude that both J and eh{-^w) admit minimiz- 
ers, their uniqueness is equivalent to the uniqueness of the solutions of the 
corresponding Euler-Lagrange equation ([3]). If uniqueness of the solution is 
satisfied, then the algorithm ([s]) can be reformulated equivalently as the fol- 
lowing two-step iterative procedure: 

• Find u'^'^^^\ which solves 

(y;(n)(^(n+l))V + M£(^("+i) _ g)v)dx = Vt; G W^'^U' R); (6) 
A 

Compute directly ^i;("+^) by 






w^-^'^ = ehV ,, l...,, A-:=l ' 6, ' H jj^j;^ < eh < }^ 



|(u(«+i))'r eh ■ I K"*"^'')' 



eft. 



otherwise. 



In the following sections we illustrate the finite element approximation of the 
Euler-Lagrange equation ([s]) similar to [151 section 8]. However, the inter- 
esting solutions may be discontinuous. In order to deal promptly with the 
discontinuity, we have studied the apphcations to gray-scale 1-dimensional 
inpainting problem by discontinuous Galerkin method for total variation min- 
imization. 

3. Finite element method for total variation minimization 

3.1. Finite element method formulation for problem M). 

Denote A = -^y^, then for a given gradient weight w^'^\ the finite element 
method for solving (pj) is to find m("+^) such that 

a(M("+i), v) =<F,v> ^ve W^-\Q] R), (7) 



where 



Jq 



and 

< F,v >= / Xgvdx. 
Jq 

Suppose the problem domain Q is discretized into A^ equal size of elements: 

= xo = xi < ■ ■ ■ < xn = i, denote Im = {xm, Xm+i) and h the mesh size. 

The integral for the m^^ element is 

"Xm + l ^ 

The trial function u is expressed as 
with the usual nodal basis functions 



'?^y^) = — : — ' M^) 



h ' ^^' ' h 

In our example (Section |5|, the value of g in each element is a constant (we 
denote it by g) and the value of A is either or ^. Now we could compute 
the element matrix and the element load vector 

Assembling the element matrices and element load vectors, we could obtain 
the linear system A^^+^^u^^+i) = b("+^). 
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3.2. Numerical implementation of the alternating-minimization 
algorithm. 

Input: Data vector g, eh > 0, initial gradient weight w^'^^ with eh < 
^(0) < J_^ number rimax of outer iterations. 

^h ^ 

Parameters: Positive weight A. 

Output: Approximation u* of the minimizer of Q. 

u(0) := 0; 

for ra := to rimax do 

Compute u("+^) such that A^^+i^u^'^+i) = b; 

Compute the gradient (u("+i)|,„)' = m^"+'Vi' + ^^S^Ii V2' = -^ + ^; 

" I («(" + !))' I fh' 

endfor 



M* := m("+^). 

4. Discontinuous Galerkin method for total variation minimization 

4.1. Discontinuous Galerkin method 

In this section, we will use Discontinuous Galerkin method to solve the 
same problem computing the solution of Euler-Lagrange equation (pi). Let 
us consider the problem 

- {wu'Y + \u = \g (8) 

Let = xo = xi < • • • < xat = 1 be an uniform partition, denote In = 
{xn,Xn+i)- Denote by Vi the space of piecewise discontinuous polynomials 
of degree 1: 

^i = {v: v\j„ e Pi(/„) Vn = 0, . . . , AT - 1}, 

where Pi{In) is the space of polynomials of degree 1 on the interval /„. 
Then we define the jump and the average of v at the boundary points of 

In. 

[f (x„)] = v{x-) - t;(x+), {v{xn)} = -{v{x-) + v{x^)) \fn = l,...,N-l. 
We also extend the definition of jump and average at Xq and x^: 

[v{xo)] = ~v{xq), {v{xo)} = v{xq), [v{xn)] = v{x~^), {v{xn)} = v{xj^). 



Next we introduce the penalty terms of the solution: 

N 
Jo{u,v) = '^T[u{Xn)][v{Xn)]. 



n=0 

where a is the real nonnegative number and h is the mesh size. 

Now we multiply ([s]) by f G Pi and use integrating by parts on each 
interval /„: 

Xn + l ^ _ /-Xn + l ^ 

{wu'v' + Xuv)dx — wu'v\ ++^ = / Xgv. 
By adding all N equations above, we obtain 

y I {wu'v' + Xuv)dx — y^ wu'v\ X^^ = / Xgv. 

n=0 "^^^ n=0 " •^'^ 

Then we have 

y / {wu'v' + Xuv)dx — \2{w{Xn)u'{Xn)}[v{Xn)]= / XqV. 

n=0 "^^" n=0 "'0 



If M is a solution of (pi), then u is continuous( [^(xn)] = for all 1 < n < 
A^ — 1 ), thus u satisfies 

N-l „a. , , N N 

y^ / (wmV + Xuv)dx - y^^{w{Xn)u' {Xn)}[v{Xn)\ + l3^Y^{w{Xn)v' {Xn)}[u{Xn)] + Jo{u,v) 

n=0 "^^^ n=0 n=0 

Xgv + /3(-u;(a;o)t^'(a;o)u(xo) + w(xAr)t;'(xAr)M(xAf)) + T-(M(xo)f (xq) + M(xAr)t;(xAr)). 
h 



Now the DG methods for solving ([8j) is to find u E Vi such that 

a(u,'y) =<F,w > V-y Gl^i, (9) 

where 



N-l ,,,,. , , N 



a(M,t>)=N^ / {wu'v' + Xuv)dx — y \w{Xn)u'{Xn)}[v{Xn)\ 

n=0 "^^'1 n=0 



n=0 



is the DG bilinear form, and 

f^~ a 

<F,- >= / \gv+l3(-w(xo)v\xo)u(xo)+w(xN)v'(xN)u(xN)) + Ti'^i^o)'i^ixo)+uixN)vixN)) 
Jo h 



is the hnear form. In our example (Section pi), we take the parameter (3 = 1 
so that the DG biliear form is symmetric. 

4-2. Linear system 

In this subsection, we derive the linear system obtained from the DG 
method. We choose for local basis functions of Pi{In) the nodal basis func- 
tions, i.e, Pi(/„) = span{0",02} with 



The global basis functions {$"} for the space Vi are obtained from the local 
basis functions by extending them by zero: 

1 otherwise. 
We can then expand the DG solution as 

N-l 2 
m=0 j=l 

Inserting this form of u^'^ into the scheme d9]), we get 

N-l 2 

^^afa($f,$n=<F,$^>, VO<n< Ar-l,Vl <2<2. 

m=0 j=l 

We would then obtain a linear system Aa = b, where a is the vector with 
components aj^, A is the matrix with entries a($'^, ^^), and b is the vector 
with the components < F, $" >. 



4-2.1. Computing the matrix A 

In this section, we will first show how to compute the local matrices. 
We will regroup the terms a(<l>™, $") into three groups: the terms involving 
integrals over J„, the terms involving the interior nodes x„, and the terms 
involving the boundary nodes Xq and x^. 

Firstly, we consider the term corresponding to the integrals over /„. On 
each element J„, the DG solution u^'^ can be expressed as 

u^^{x) = <0^(x) + a^(j)^{x) Vx e /„. (11) 



^n + 1 , , /"^n + l 



Thus, using (|11 ) and choosing v = cf)^ for i = 1,2, we get 

j—l J Xn 

This linear system can be written as A„a;", where 



wiu^"") (0r)'+Aw^^C)rf^ = y.<^l I {w{r,)'{€)'+H';€)dx W = 1, 2. 






n 



a 



2 



We could compute the A„: 

^"~ /i V -1 1 y* V 1/6 1/3 

Second, we consider the terms involving the interior nodes a;„. Let us 
express 

- {w{Xn){u'''')\Xn)}HXn)]+P{w{Xn)v\Xn)}[u''''{Xn)] + ^ [w^^(x„)] [t;(x„)] 
^n ~r Ctt, ~r (In \ Cn? 

where the terms are defined as below: 

hn = \n^{xt^{u^^)\xtM<) - ^wix^^^^ixtWixt) + ^u^^ {x^W (x^) , 

Cn = --W{x-){u^^y{x-)v{x-) + -W{X-)U^^{X-)V'{X-) + -U^^ {x;,)v' {x') , 

dn = -^^(x+)(«^^)'(x+)t;(x;) - ^w{x-)u^^{x^^)v\x-) - ^u^^{xt)v'{x-), 
en = lw{x-){u^^y{x-)v{x^) + ^w{x^)u^^{x-)v'{x^) - lu^'^ix-Wix^J. 



Now with the expression (11) and the choice f = 0", the four terms defined 



above yields the local matrices B„, C„, D„ and E„ which are: 

l_ ( -w{xt) + Pw{x+) + 2a w{x+) 
2h V -(3w{x+) 

1 ( -Pw{x-) 

2h \ w{x^^) -U7(a;+) + (3w{x1^) + 2a 

2/i V ^(x+)-/3w(x„)-2a -w(a;+) 
1 f -w{x-) w{x-) - Pw{x+) - 2a 



B„ 

D„ 

^" ~ 2/1 V /3«;(a;^ 

Finally, we compute the local matrices from the boundary nodes xq and 



Xn- 



a 



/o = w{xq){u'^'^)' {xq)v{xq) - (3w{xo)u^'^{xo)v'{xo) + '^u^'^{xo)v'{xo), 

n 



JN = -w{xN){u^^y{xN)v{xN) + f3w{xN){u^'^){xN)v' (xn) + -U^^ {x n)v' {x n) , 

which yields the local matrices Fq and F^r: 

p _ -1 / -w{xq) + (5w{xq) + a w{xq) 
° ~ /i V -Pw{xo) 

_ —I5w{xn) 

N 



h \ w{xn) —w{xn) + I3w{xn) + a 



Assuming that the unknowns are listed in the following order: 



fa?, ao, a],a\, a?, a^ 






'1) "25 "1)"25 "15 "2? • • • ? "1 )"2 >) 

we obtain the global matrix A which is block tridiagonal 

/ Mo Di \ 

El M D2 



V 

where 

M = A„ + B„ + C„+i,Mo 



E^_2 M T>N-i 



Ao + Fo + Ci,M 



N 



^N~l 



' N 



c 



N-l- 



4-2.2. Computing the right hand side b 

Each component of b can be obtained by computing 

<F,<I>^>= / Xg¥^^dx + l3{-w{xo)i(^iyixo)u{xo) + w{xN)i(^nyixN)u{xN)) 
Jo 

+ -{u{xo)^i,{xo) + uixNWni^N))- 

n 
Because of the local support of ^l^, the first term is reduced to 

/ \g^\dx= j \9(l)ndx. 

Jo J Xn 

We arrange the components of b in an order consistent with the order of the 
unknowns a": 

l"H "2' "D "2) "l' "2' • • • ' "l '"2 /' 

where the first two components and last two components are 

_ Xgph f3w{xo)u{xo) au{xo) 

^'~ 2 ^ ^ h ^ h ' 

_ Xgph _ (3w{xo)u{xo) 

2 h 

,N-i _ ^9N-ih (3w{xn)u{xn) 
"1 — 



2 h 

^' - 2 ^ h + 

and the other 2(A^ — 2) components are 



^_i _ XgN-ih /3w{xn)u{xn) olu{xn) 
2 " 2 + h ^ h ■ 



ya^^9nh vi< n < A^ - 2, Vl< i < 2. 

t 2 ' — — , _ _ 

5. Numerical examples and future study 

In this section we show numerical results of the applications of the fi- 
nite element method and Discontinuous Galerkin method for total variation 
minimization described above. 

We would first use finite element method for total variation minimization 
the do some experiments with different values of parameter A in the alternat- 
ing minimization algorithm. Let us consider the same signal with the same 
inpainting interval (|, |), see Figure ^1. We conclude that quality of the 
recovery image is becoming better as the value of A increases. 

10 
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Figure [5jl The three 'pictures from left to right represent the result of three 
different values of X (X = 10, 100, lOOOj, respectively. 

Second, we fix the outer iteration to 20 and compare the convergence 
speed for three different values of A (A = 10, 100, 1000), see the left picture 
of Figurej5J2. And we conclude that the convergence speed increases as the 
value of A increases. 



Finally, we modify w^^+i) from e^V ,, , i,,,,. A— to (e/iV ,, , i ,,,,, , , 



-A 



1 \2- 



and check the convergence speed for different value of r(r=l, 0.9, 0.8, 0.7, 0.6, 0.5). 
The result is shown in the right picture of Figure |5j2. And we conclude that 
the convergence speed increases as the value of r decreases. 



Convergence speed with different luii 




Convergence speed with different tao 




Figure |5j 2^ r/ie left picture illustrates the convergence speed for three different 
values of X (X = 10, 100, lOOOj; the right picture shows the convergence speed 
for different value of r (t=1, 0.9, 0.8, 0.7, 0.6, 0.5). 

Next we consider a signal with a jump (see Figure [5 
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Figure |5|3 The signal of a step function. 



Let us recover the signal(Figure p|3) by finite element method and dis- 
continuous Galerkin method for total variation, respectively. The results are 
shown in Figure [5J4. We observe that the finite element method for total vari- 
ation minimization couldn't preserve the jump very well from our example. 
However, the discontinuous Galerkin method for total variation minimization 
preserves the jump rather well. 

Our future study aims at the construction, analysis and implementation 
of new adaptive discontinuous Galerkin (DG) solvers for total variation min- 
imization problems in two space dimensions. These methods are based on 
re-weighted least squares and are implemented by nester outer and inner it- 
erations. The adaptivity concerns not only the space discretization but also 
the parameters involved in the inner and outer iterations as well as in the DG 
discretization. The inner iteration should be robust with respect to both the 
discretization and the gradient weights, and the number of inner iterations 
should be controlled in proper way. The robustness property is obviously 
connected with the preconditioning. 



200 elements: iim=1000; 20 iterations 







- 


- 



200 elements: Mu=100: penalt}'= 20: 20iterations 
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Figure |5}4 The left picture shows the recovery image by the finite element 
method for total variation minimization; and the right picture illustrates dis- 
continuous Galerkin method for total variation minimization. 
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